Instability of defensive alliances in the predator-prey model on complex networks 
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A model of six-species food web is studied in the viewpoint of spatial interaction structures. Each 
species has two predators and two preys, and it was previously known that the defensive alliances of 
three cyclically predating species self-organize in two-dimensions. The alliance-breaking transition 
occurs as either the mutation rate is increased or interaction topology is randomized in the scheme of 
the Watts-Strogatz model. In the former case of temporal disorder, via the finite-size scaling analysis 
the transition is clearly shown to belong to the two-dimensional Ising universality class. In contrast, 
the geometric or spatial randomness for the latter case yields a discontinuous phase transition. The 
mean-field limit of the model is analytically solved and then compared with numerical results. The 
dynamic universality and the temporally periodic behaviors are also discussed. 
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The rock-scissors-paper (RSP) game pj gives a typ- 
ical three-strategy model of cyclical predator-prey food 
chain if rock, scissors, and paper are replaced by differ- 
ent competitive species in biology. In this simple game, 
a rock beats a pair of scissors, a pair of scissors beat a 
sheet of paper, and a sheet of paper beats a rock. As 
generalizations of RSP game, the food webs composed of 
different number of species show a variety of interesting 
behaviors 2]. Especially, the models proposed by Szabo 
et al. in Refs. Q and Q have been shown to self-organize 
to form the defensive alliances, which becomes unstable 
as the mutation rate is increased. 

In the present study, we use the model A in Ref. |3( 
described by the interaction topology in Fig.^a), where 
the six species weave a food web that is more complicated 
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FIG. 1: (a) The interaction topology used in this work (Model 
A in Ref. 0). Arrow e.g., from to 1 indicates that the 
species is a predator to the species 1 (and thus the species 
1 is the prey to the species 0). As time goes on the ini- 
tial random configuration in (b) evolves to (c), where only 
three species [(0,2,4) or (1,3,5)] exist to form the defensive 
alliance. Different colors indicate different species in (b) and 
(c). (Color online.) 
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than the RSP game: Each species has two predators, 
two preys and one neutral interacting partner. There 
indeed exists ecological systems showing the cyclic domi- 
nance: Examples include the three-morph mating system 
in the side-blotched lizard in which each morph domi- 
nates another morph when rare || . Another example is 
the system of three different populations of Escherichia 
coli composed of toxin-producing (T), toxin-resistant 
(R), and toxin-sensitive (S) strains. If the growth rate 
of each group satisfies S > R > T, the RSP game cap- 
tures the cyclic dominance in the system: T dominates 
S (S is killed by T), S dominates R due to the higher 
growth rate, and R dominates T from the same reason. 
If bacteria produce two different toxins, the food web is 
constructed from nine different species as has been well 
studied in Ref. 0. 

The main interest here is to study the effect of the spa- 
tial interaction structure and for that purpose we play 
the predator-prey game on the complex network struc- 
ture of the Watts-Strogatz (WS) network [6j constructed 
on two-dimensions (2D) as follows: (1) We first build 
the 2D L x L (N = L 2 ) regular square lattice. (2) Ev- 
ery bond is visited once and with the rewiring probabil- 
ity a is rewired to the randomly chosen other site. The 
above procedure then yields a network structure which 
possesses characteristics such as the short characteristic 
path length 0,0- 

Once the network is built as described above, the time 
evolution of the system obeys the following rules Q: (1) 
Six species are scattered randomly on a square lattice as 
in Fig.^b), then (2) the species on each randomly cho- 
sen site is mutated to one of its predating species with 
the mutation rate P. (3) If no mutation occurs (with 
the probability 1 — P), one of the nearest neighbors is 
chosen to interact, and the dominant one survives and 
invades the subordinated one. For example, if the pair 
of species and 1 are chosen, the species 1 is replaced 
by the species [see Fig. IHa)]. If two neutral partners 
has been chosen, i.e., no arrow connects the two species 
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in Fig. ma), nothing happens. For the 2D regular square 
lattice corresponding to the rewiring probability a = 0, it 
has been found that the defensive alliances composed of 
three species, (0,2,4) or (1,3,5), are spontaneously formed 
at a small mutation rate and the other species that does 
not belong to the alliance dies out Q [see Fig.EJc)]- As 
the mutation rate is increased, it has been shown [3] that 
the defensive alliance becomes unstable and there occurs 
well-defined phase transition of the universality class of 
the 2D Ising model @ , which is expected because of the 
existence of the Z 2 symmetry: Interchange of two al- 
liances (0,2,4) <-> (1,3,5) do not change the game rules, 
which we call here inter-alliance symmetry. The sponta- 
neous breaking of the defensive alliance originates either 
from the high mutation rate P or from the high degree of 
structural randomness controlled by the rewiring proba- 
bility a. In other words, the instability of the defensive 
alliance is induced either by the temporal randomness or 
by the structural randomness. In reality, these two differ- 
ent types of randomness may coexist. However, from the 
practical computational difficulty we in this paper only 
investigate the effect of each randomness separately. 

We first investigate the phase transition for the 2D 
regular square lattice corresponding to the WS network 
with the rewiring probability a = 0. Although this was 
previously studied 0, we in this work use the extensive 
finite-size scaling analysis to confirm not only the 2D 
Ising static universality class but also to identify the dy- 
namic universality class detected by the dynamic critical 
exponent z. In order to describe the alliance breaking 
transition we define the order parameter which functions 
as a magnetization in the Ising model: 



m = (c + c 2 + c 4 ) - (ci + c 3 + c 5 ), 



(1) 



where c s = N s /N (s = 0, • • • , 5) is the density of species s 
with N s being the number of sites occupied by s. When 
P is sufficiently small, the defensive alliance is formed 
with Co ~ C2 ~ C4 Ri 1/3 and c\ ~ C3 w C5 « or 
vice versa, leading to m ~ ±1 (ordered phase). When 
the system recovers its full symmetry at a high mutation 
rate, all species have the same density, yielding m « 0. 
Only for convenience we use 



M = ln(l/P) 



(2) 



as a control parameter. High mutation rates correspond 
to the small values of the mutation parameter fx, and 
thus qualitatively speaking, the physical meaning of ll 
resembles that of the inverse temperature in the standard 
statistical mechanics. 

Numerical simulations are performed on systems of 
sizes varying from L = 32 to 192 with the periodic 
boundary conditions employed. After equilibration for 
10 4 steps per site, which is sufficiently long enough, the 
thermal average is computed for later 10 4 steps at least. 
The measured quantities are (|m|), the Binder's cumu- 
lant @ 



and the susceptibility 

x = N({m 2 )-(\m\) 2 ), 



(4) 



where (...) denotes the thermal average. In order to study 
dynamic critical behavior, we also measure the autocor- 
relation function as a function of time t defined by 



C{t) = (\m{t)m(0)\) - (\m\) 2 



(5) 



U L = 1 



(m 4 ) 
3(m 2 ) 2 



(3) 



The finite-size scaling forms of measured quantities are 
written as 

(\m\)=L-e/ v m((Li- f i c )L 1 /»), (6) 
U l = U((li-ll c )L 1 ^) 1 (7) 

x = l^x((m-Mc)£ 1/,/ ), (8) 

C(t)/C(0) = C{tL- z ), (9) 

where m, U, x, and C are suitable scaling functions, ll c is 
the critical value of li, and /?, 7, v are standard critical 
exponents |8j while z is the dynamic critical exponent 
defined at li c from the divergence of the relaxation time 
scale (r ~ L z ). 

Figure [21 summarizes the numerical results for the 
phase transition in the 2D regular lattice. The order pa- 
rameter (\m\) shown in Fig. Efa) exhibits the existence 
of the transition, which is analyzed in detail in Fig. |2Ib) 
and (c) by using the finite-size scaling form in Eq. JJJ. 
The critical point ll c = 6.50(4) as well as critical expo- 
nents (3 « 1/8 and v w 1 are obtained, which are con- 
firmed again from the finite-size scaling of the Binder's 
cumulant shown in Fig.[3fd). The divergence of the sus- 
ceptibility in Fig. Efe) is analyzed to get 7 w 7/4. At 
H = Li c , we compute the autocorrelation function Q and 
plot it in Fig. ffif) in accord with the scaling form in 
Eq. 10: All curves at different sizes collapse well to a 
single curve with the dynamic critical exponent z ss 2. 
All these findings clearly confirm that the alliance break- 
ing transition induced by the mutation belongs to the 2D 
Ising universality class as was known in Ref. 

We next study the phase transition induced by the spa- 
tial randomness introduced via the long-range shortcut 
in the WS network model. One can motivate the study 
along this direction since in real systems, the spatial 
interaction topology among species can be much more 
complicated than the nearest-neighbor interaction on a 
regular square lattice. Distinguished from the regular 
network, the small- world network 0, Q , which captures 
characteristics of many real networks very well, has a re- 
markably small average path length similar to the glob- 
ally coupled network (or the mean-field case). The three- 
strategy RSP game has been studied on some small- world 
networks and the periodic flourishes of three strategies 
were found to happen during the time evolution |llj| . 

We fix the mutation parameter to Ll = 7.0, which is 
well inside the ordered phase with the defensive alliance 
formed for the 2D regular square lattice (see Fig.[2J). Use 
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FIG. 2: Phase transition in the 2D regular lattice in terms 
of the mutation parameter /i. (a) The order parameter (|m|) 
versus /i clearly shows the existence of the alliance breaking 
transition, (b) A unique crossing at fi c — 6.50(4) with /3/v « 
1/8 is obtained from the finite-size scaling form of (|m|). (c) 
All data points for (\m\) in (b) collapse into to a smooth 
curve by using scaled variables, (d) Binder's cumulant at 
different sizes cross at \x c = 6.50(4) in agreement with (b). 
The inset in (d) shows that AU L ~ L 1/v to yield v « 1, 
where AUl = Ul(j-Ii) — Ul{^2) with /ii^ are two adjacent 
points near /j, c . (e) The susceptibility \. The inset shows 
a log- log plot of x versus L at fi c leading to 7 ~ 7/4. (f) 
The autocorrelation function C(t)/C(0) at fi c versus tL~ 19 
with z — 1.9. All curves for different sizes collapse well to 
a single curve, indicating that the dynamic critical exponent 
z = 1.9(2). (Color online) 




FIG. 3: Phase transition in the WS network in terms of the 
rewiring probability a. The mutation parameter is set to /1 = 
7.0. (a) The order parameter (\m\) as a function of a shows 
a sudden drop at a c « 0.28. The change of (\m\) becomes 
steeper as L is increased, indicating that the transition is 
discontinuous one. (b) The normalized histogram H(m) of m 
for L = 1024 is displayed at a = 0.27, 0.28, and 0.29. The 
abrupt change in the form of the histogram between a — 
0.28 and 0.29 clearly confirms again discontinuous transition. 
(Color online) 



FIG. 4: Time evolution of density of species c 3 (t) for the 
WS network for L = 1024 at fj, = 7.0. (a) At a = 0.27, 
the inter- alliance symmetry is broken, suggesting tn/0, and 
each species in the alliance dominates all the other species 
in a time-periodic fashion. It is interesting to note that the 
non-dominant alliance is also formed and even the member of 
the non-dominant alliance prevails periodically, (b) At a = 
0.29, the system recovers its full symmetry and each species is 
equivalent to others. However, the time-evolution of c s {t) is 
still periodic, in contrast to the disordered phase in 2D regular 
lattice. The time t is measured after stationarity is achieved. 
(Color online) 



of a very large value of /x, corresponding to very small mu- 
tation rate, was found to make the system approaches a 
local dynamic fixed point and then the system stays there 
forever. At /i = 7.0, which is small enough to ensure the 
equilibration and large enough to make the system or- 
dered at small a, the system is found to undergo a phase 
transition at a = a c « 0.28(1) as displayed in Fig. [3] 
The abrupt drop down of the order parameter (|m|) at 
the transition as displayed in Fig. OJa), together with 
the change of the histogram Him), normalized to sat- 
isfy ^ m H{m) = 1, in Fig. 01b), clearly indicates that 
the transition is of the discontinuous nature in a sharp 
contrast to the finding for 2D regular lattice at a — 
(see Fig. 0). In more details, the sudden change of the 
peak position of H(m) between a — 0.28 and 0.29 in 
Fig. EJb) is interpreted as a strong evidence of a dis- 
continuous transition. Continuous transition in general 
exhibits the continuous shift of the peak position toward 
m = as the critical point is approached from the or- 
dered phase. 

In general, one can study the phase diagram of the 
model in the 2D parameter space of (/x, a). Due to the 
practical difficulty to cover the whole parameter space, 
we in this work only explore the phase transitions along 
the two straights line in the (/x, a) plane: One on the 
axis (/i, a — 0), and the other on the line (/i = 7.0, a). 
Both parameters \x and a control the amounts of ran- 
domness (temporal and spatial ones, respectively), and 
accordingly the ordered phase with the defensive alliance 
becomes unstable as either [i is decreased (i.e., the mu- 
tation rate P is increased) or a is increased. Not only 
the phase transitions have different natures (continuous 
one belonging to the 2D Ising universality class for the 
former, and discontinuous one for the latter), but also 
the ordered and disordered phases in each case are very 
much different in terms of the time evolution of densi- 



4 



0.8 
0.6 
0.4 
0.2 








f 

1 j 


N= 32^ ' 




64^ 




128; - 




256; " 




512 



0.1 



(b) 



10 15 20 



FIG. 5: Numerical results for the mean-field case, (a) The 
order parameter (|m|) versus \i for N = 32 2 , 64 2 , ■ ■ ■ , 512 2 . 
As N is increased the region for the ordered phase shifts to- 
wards higher values of fi, indicating that the defensive alliance 
vanishes in the thermodynamic limit at any nonzero value of 
the mutation rate, (b) Time evolution of the order parame- 
ter m(f) is not periodic at all, different from the WS network 
shown in Fig. 2] (Color online) 



ties of species. In the ordered and disordered phase on 
the axis of a = 0, the density of each species does not 
fluctuate much but stays at almost the same level: For 
ordered phase at fi > ^i c , c s w 1/3 for s 6 {alliance} 
and c s w otherwise, while for the disordered phase 
at \i < /i c , c s ~ 1/6 for all species. In contrast, the 
time evolution c s (t) for the case of the WS network at 
(/i = 7.0, a) is strikingly different. In Fig.^Ja), the time 
evolutions of densities of species are shown for /i = 7.0 
and a = 0.27(< a c ). It is clearly shown that the inter- 
alliance Zi symmetry is broken, indicating that the sys- 
tem has a nonzero value of the order parameter. Another 
very important observation one can make from Fig. Ufa) 
is that each species within the alliance (s = 0, 2, 4) cycli- 
cally dominates all the others in a very regular way, and 
very interestingly, the species that does not belong to the 
dominant alliance also prevails in a time-periodic man- 
ner. Even in the disordered phase at a > a c , the periodic 
dominance persists while the inter-alliance symmetry is 
fully recovered [see Fig. Efb)]. 

We finally investigate the mean-field limit of the 
game where all individual species interact with all 
the other species in the system. The master equation for 
the number N s of species s is given by 



AN, = P - 



2 2 

+ 2(1 - P)c s ( Cll +c t 



^71 



cj a ), (10) 



where AN S = N s (t + 1) — N s (t), i\ and i^ (ji and J2) 
are two preys (two predators) of s. For example, the 
species s = has i\ — 1, ia = 2,ji = 4,j2 = 5 [see 
Fig. Ola)]. The term proportional to P describes the 
decrease of N s by the mutation from s to its predators 
and increase of N s by the mutation from its preys to s. 
The other term originates from the interaction of s to 
other species: When s meets its preys (predators) N s is 
increased (decreased). The numerical factor two in front 



of (1 — P) is due to the interaction of other species with 
s. If we start from the situation when cq = c 2 = C4 = 
(1 + m)/6 and c\ = C3 = C5 = (1 — m)/6, with the order 
parameter m in Eq. |JI|I. the above master equation is 
reduced to a very simple form 



dm 

= m(t + 1) — m(t) 



-NPm, 



yielding the solution 

m(f) = m(0) exp(— f/r) 



(11) 



(12) 



with r = PN. The mean-field solution indicates that 
there is no ordered phase at any nonzero value of the mu- 
tation rate, i.e., the defensive alliance cannot be formed in 
the mean-field limit. It is also noteworthy that the time- 
periodic behavior observed for the WS network in Fig. 0] 
ceases to exist in the mean-field case. The simulation re- 
sults displayed in Fig. for the mean-field model are in 
perfect agreements with the above analytic findings: The 
ordered phase that appears to exist for small system sizes 
drifts away toward the region of higher value of /1 as A^ 
is increased, suggesting the disappearance of the ordered 
phase in the thermodynamic limit [see Fig. Efa)]. There 
is no time-periodic behavior of m in equilibrium as shown 
in Fig.|5fb), where f is measured after equilibration. 

In summery, we have investigated the instability of the 
defensive alliances for the simple food web of six species 
in three different spatial interaction structures: the 2D 
local regular square lattice, the WS network, and the 
globally coupled mean-field network. The temporal ran- 
domness imposed by the mutations as well as the spa- 
tial randomness in the interaction structure tuned by the 
rewiring probability in the WS network has been shown 
to make the defensive alliance unstable. When the muta- 
tion rate is increased for the 2D square lattice the alliance 
breaking transition has been clearly identified to belong 
to the 2D Ising universality class due to the common Z2 
symmetry. On the other hand, when the rewiring proba- 
bility is increased, the transition becomes discontinuous, 
and around the transition the natures of the ordered and 
disordered phases are very different from the 2D square 
lattice. The mean-field model has also been studied ana- 
lytically and numerically with the results that the defen- 
sive alliance cannot be formed at any value of the mu- 
tation rate and that the time-periodic behavior observed 
in the WS network is not seen any longer. 
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